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CHAPTER  I 


INTRODUCTION . 

Many  types  of  structures  are  designed  for  use  below  and/or 
above  the  surface  or  on  the  shoreline  of  oceans  and  various  large 
inland  water  areas.  They  can  be  floating  and/or  structurally  supported 
to  withstand  gravitational  forces.  These  structures  are  constructed  to 
withstand  wave  forces  to  the  extent  present  day  knowledge  of  water 
wave  mechanics,  standard  of  practice,  and  appropriate  factor  of  safety 
allow.  A criterion  used  in  design  may  be  the  "design  wave",  a 
statistical  culmination  of  known  wave  heights  recorded  in  the  area 
combined  with  an  appropriate  crest  length  which  is  used  to  compute  the 
largest  horizontal  forces  which  the  structure  will  be  required  to 
withstand. 

An  important  design  assumption  used  for  many  structures  built 
to  withstand  water  wave  forces  is  the  reflection  coefficient.  Defined 
as  the  ratio  of  the  height  of  the  wave  reflected  from  the  structure  to 
the  height  of  the  incident  wave  coming  in  contact  with  the  structure 
{!}*,  it  is  an  indicator  of  wave  energy  reflected  from  the  structure 
as  the  square  root  of  wave  energy  is  directly  proportional  to  wave 
height  or  amplitude.  This  reflected  energy  is  an  indicator  of  the 
direct  effect  the  structure  will  have  on  any  neighboring  structure  or 
shoreline.  Combined  with  knowledge  of  the  transmitted  energy  which 
passes  through  the  structure  and  possible  motion  of  the  structure 
itself,  it  can  be  used  to  determine  energy  absorbed  by  the  structure 
and  forces  which  the  structure  will  be  required  to  withstand  to 
perform  properly. 

^Numbers  shown  in  brackets  refer  to  the  References  listed  on  pp.  39-^0. 


2 


Specifically,  breakwaters,  floating  bridges,  piers,  wharves, 
and  other  structures  designed  to  withstand  and  even  modify  water 
wave  forces  have  a particular  requirement  for  knowledge  of  reflection 
coefficients  both  in  their  design  and  for  xollow-up  monitoring  of  true 
effects  after  construction  is  completed.  Determination  of  reflection 
coefficients  either  from  model  study  during  design  or  from  study  of 
the  finished  structure  itself  is  not  always  simple,  even  if  the 
situation  allows  measurement  of  both  the  incident  wave  sea  state 
alone  and  the  combined  reflected  and  incident  wave  sea  state.  These 
measurements  are  highly  dependent  on  angle  of  attack  of  the  incident 
wave  train,  local  conditions  (water  depth,  etc.),  and  the  particular 
frequencies  and  amplitudes  associated  with  different  incident  wave 
trains  found  at  different  times.  The  problem  is  particularly 
difficult  when  one  is  working  with  a structure  such  as  a floating 
bridge  where  it  is  impossible  to  measure  the  incident  wave  sea  state 
prior  to  its  combination  with  reflected  wave  components. 

Morden  {2,  3},  Thornton  and  Calhoun  {4}  and  others  have 
developed  deterministic  approaches  to  computing  reflection  co- 
efficients from  the  combined  incident/reflected  sea  state  based  on 
linear  deep  water  wave  theory.  Their  methods  require  the  use  of  two 
or  more  wave  measurement  devices  in  the  combined  wave  field. 

In  contrast,  the  thrust  of  the  current  investigation  is 
the  separation  of  incident  and  reflected  wave  energy  _rom  the  combined 
sea  state  using  only  one  wave  measurement  device  and  cepstral  analysis 
techniques.  The  present  use  of  cepstral  analysis  is  primarily  in 
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seismic  epicenter  estimation,  speech  parameter  estimation,  photo- 
graphic image  enhancement,  radar,  sonar,  and  other  fields  where  one 
desires  separation  or  identification  of  the  reflected  signal  from  the 
original  signal  or  the  time  of  delay  between  the  two  signals.  The 
literature  contains  several  references  to  recent  capstral  analysis 
work  performed  in  the  speech  and  seismic  arenas  but  apparently  no  one 
has  ever  attempted  application  of  this  process  to  water  waves. 

Cepstral  analysis  wave  recovery  and  decomposition  techniques 
assume  neither  the  reflection  coefficient  nor  the  delay  time  are 
functions  of  frequency.  As  this  is  not  an  absolutely  correct 
assumption  with  water  waves,  it  is  probably  understandable  why  this 
technique  has  not  been  actually  attempted  with  water  waves.  However, 
cepstral  analysis  decomposition  would  allow  determination  of 
reflection  coefficients  using  only  one  wave  measurement  device 
collecting  data  from  the  combined  incident/reflected  wave  field.  It, 
therefore,  would  be  of  great  benefit  if  this  method  could  be  reliably 
applied  with  respect  to  ocean  waves. 

Chapter  Two  discusses  very  briefly  some  of  the  general 
mathematical  concepts  employed  in  this  work  and  directs  the  reader  to 
other  references  if  one  whishes  to  delve  into  discrete  time  series 
analysis  in  more  detail.  It  by  no  means  is  an  exhaustive  explanation 
of  background  material  but  is  instead  a hint  at  what  the  field 
involves. 

Chapter  Three  contains  a description  of  cepstral  analysis 
in  general.  Starting  with  a description  of  the  power  cepstrum,  it 
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then  discusses  the  complex  cepstrum  and  the  process  of  incident  wave 
recovery.  The  chapter  closes  with  limitations  anticipated  when 
applying  cepstral  analysis  to  ocean  waves. 


Chapters  Four  and  Five  discuss  the  test  models  used  in 
applying  cepstral  analysis  to  simulated  water  waves  and  the  results 


of  this  application.  The  first  model  is  a finite  summation  of  sinusoid 
signals  containing  random  phase  components  keeping  both  delay  time 


and  reflection  coefficient  constant  for  all  frequencies.  The  second 
model  is  the  same  summation  of  sinusoid  signals  as  the  first  but  the 
delay  time  of  the  reflected  wave  is  allowed  to  vary  with  frequency. 
The  third  model  is  a summation  of  sinusoid  signals  allowing  the 
reflection  coefficient  as  well  as  the  delay  time  of  the  reflected 
wave  to  vary  with  frequency. 

Chapter  Six  is  a brief  look  at  some  real  world  ocean  data 
using  cepstral  analysis.  A hint  is  given  of  problems  remaining  to 
be  solved  before  this  technique  can  be  easily  applied  in  the  field. 

Chapter  Seven  lists  the  conclusions  from  this  work 
including  recommendations  for  further  study.  Overall,  the  process 
shows  promise  and  with  more  research  could  become  a valuable  tool 
in  the  decomposition  of  combined  incident/reflected  sea  states. 
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CHAPTER  II  . 

BACKGROUND  INFORMATION 

The  work  throughout  this  study  is  based  on  use  of  the 
discrete  Fourier  Transform,  specifically  the  FFT  or  Fast  Fourier 
Transform.  Reference  will  be  made  for  illustrative  purposes  to 
other  transform  methods  but  all  analysis  was  actually  performed  using 
the  FFT.  For  a general  review  of  the  FFT  and  what  it  means  and  can 
do,  the  reader  is  referred  to  Bergland’s  article  "A  Guided  Tour  of 
the  Fast  Fourier  Transform"  {5}. 

Although  for  simplicity  some  examples  may  be  written  in 
terminology  more  akin  to  continuous  transform  analysis,  all  work  in 
this  study  is  specifically  digital  time  series  analysis  and  is 
sometimes  only  possible  in  a sampled  time  series  sense.  Since 
Digital  Signal  Processing  or  Digital  Time  Series  Analysis  is  an 
independent  field  of  study  usually  couched  in  the  realm  of  Electrical 
Engineering,  strict  adherence  to  digital  analysis  terminology  would 
possibly  thwart  the  main  thrust  of  this  paper  which  is  to  explain  the 
cepstrai  analysis  process  In  terms  an  engineer  or  scientist  who  is 
not  thoroughly  familiar  with  time  series  analysis  may  still  understand. 
For  a detailed  description  of  the  Discrete  Fourier  Transform,  the 
reader  is  referred  to  Chapter  9 of  Stanley  {6},  Chapters  3 and  6 of 
Oppenheim  and  Schafer  {7},  Chapters  6-9  of  Brigham  {8},  and/or 
Chapter  5 of  Childers  and  Durling  {9}. 

All  computer  modeling  and  testing  was  performed  using  CSAP 
(Cepstrai  and  Spectral  Analysis  Program)  which  was  developed  by  the 
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author  from  TSAPC,  an  instructional  program  developed  by  Christensen 
at  the  University  of  Washington  used  in  teaching  spectral  analysis. 
Some  subroutines  for  CSAP  were  obtained  from  NSAP  (Numerical  Spectral 
Analysis  Program)  also  developed  by  Christensen  and  were  modified  by 
the  author  as  appropriate  to  meet  the  needs  of  this  particular 
research.  The  FFT  subroutine  used  wa3  originally  called  NLOGN 
developed  by  Thompson  and  is  described  by  Paniker  {10}.  Spectral 
computation  is  performed  using  this  FFT  subroutine  and  decimation  by 
frequency.  The  cepstral  computations  were  performed  using  a sub- 
routine developed  by  the  author  based  heavily  on  work  presented  by 
Skinner  {11}. 

Most  signal  processing  people  address  cepstral  analysis  in 
the  literature  in  terms  of  the  Z-transform.  In  this  study,  use  is 
made  of  the  Fourier  Transform  which  is  exactly  the  same  operation  as 
computing  the  Z-transfo.m  on  the  unit  circle.  In  brief,  the  Z- 
transform  is  merely  a discrete  version  of  the  Laplace  transform. 

The  Laplace  transform  (one-sided)  can  be  expressed  by 

X(s)  * /“x(t)e~stdt 

whose  inverse  is 

x(t)  “ -i— / X(s)estds 
2iri  c 

It  can  be  seen  that  if  s ■ iu  where  u =*  2:rf  the  above  relationships 
become  the  familiar  continuous  Fourier  transform  pair. 

The  Z-transform  is 

X(z)  = Z x(n)z-n 
n=o 


where  z 
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If  s ■ iiu,  then  the  Z-transform  is  computed  on  the  unit 

circle  which  yields  the  discrete  Fourier  transform, 

X(f)  - " x(n)e~iairiAt 

n»o 

or  more  correctly  as  the  finite  sum 

-i2irhn 

X(h)  - ”£*(n)e  N 
n=*o 

(where  N ■ total  number  of  sampled  points  and  h represents  digitized 
frequency  in  the  same  way  that  n represents  digitized  time) 


whose  inverse  is 


x(n) 


i NS1  X(h) 
N n=o 


i2Trhn 
e N 


In  this  way,  data  points  representing  digital  or  discrete  time  series 


operated  on  by  the  Fourier  transform  become  data  points  representing 


a digital  frequency  series,  and  the  digital  frequency  series  operated 
on  by  the  inverse  Fourier  transform  becomes  a digital  time  series. 

For  our  purposes  this  relationship  between  the  Z-transform 
and  the  Fourier  transform  is  sufficient.  The  Z-transform  is 
extensively  used  in  the  literature  as  it  is  a more  general 
mathematical  principle  that  includes  the  Fourier  transform  as  only 
a specific  part  of  its  domain.  For  a detailed  background  on  the 
Z-transform  the  reader  is  referred  to  Chapter  4 of  Stanley  {6}, 


Chapter  2 of  Oppenheim  and  Schafer  {7},  Chapter  2 of  Childers  and 


Durling  {9},  and/or  Chapter  2 of  Gold  and  Rader  {12}. 
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CHAPTER  III  . 

DESCRIPTION  OF  CEPSTRAL  ANALYSIS 

A.  Power  Cepstrum  and  Determination  of  Delay  Times. 

In  1963,  Bogert,  Healy,  and  Tukey  published  their  classic 
article  {13}  on  power  cepstral  analysis  which  describes  the  cepstrum 
(pronounced  KSp-striSm,  which  violates  a few  rules  of  English,  but 
sounds  better  than  Sep^-struin)  as  a tool  for  determining  delay  times 
between  reflected  wave  components  and  incident  wave  components 
when  recording  both  simultaneously  as  a combined  wave  field.  The 
simplest  model  for  the  combined  time  series,  y(t),  containing  only 
one  reflected  wave  component  or  echo  is 

y(t)  = x(t)+crx(t-r) 

where  x(t)  represents  the  original  or  Incident  wave  and  t represents 
the  delay  time  between  the  beginning  of  the  incident  wave  to  the 
beginning  of  the  reflected  wave.  The  a term  in  the  above  expression 
is  the  reflection  coefficient,  which  is  the  ratio  of  the  reflected 
wave  height  to  the  incident  wave  height. 

The  Fourier  transform  of  the  reflected  echo  is  obtained 
from  the  Fourier  transform  of  the  original  wave  time  series  by 
multiplying  by  ae-^1^  where  f is  in  frequency  dimensions  (Hz) . 
Therefore  the  Fourier  transform  of  the  combined  wave  time  series  is 

■ 
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Y (f ) = X(f){l+ae-2lTifT} 

Since  the  magnitude  squared  of  the  Fourier  transform  represents  the 


power  spectrum,  we  have 


| YCf  ) | 2 - |x(f)|2{l'Kx2+2acos2Trfr} 

which  represents  the  power  spectrum  of  the  combined  wave  time  series. 


Taking  the  log  of  this  power  spectrum  yields 


log|Y(f)  |2  = log|x(f)  |2+log{l+a2+2acos2irf-r} 

log|Y(f)|2  = log(l+a2)  |x(f)  |2+log(l-h-^  cos2irf  x) 

1+a2 


Using  the  log  series  expansion: 


2 3 

log(l-hx)  =*  x^L_+2L_-  •••*,|x|<l 
2 3 

we  can  expand  the  second  term  into  a convergent  series  thus: 

log(l+- 2Ct  cos27rfr)  =*  E -lull  \ .?.a-  cos27rfx)k 
1+a2  k=l  k 

The  log  of  the  power  spectrum  of  the  combined  wave  time 
series  then  contains  a cosinusoidal  ripple  associated  with  the 
reflected  wave  component.  The  "frequency"  of  this  ripple  can  be 
found  by  again  calculating  a power  spectrum.  This  second  power 
spectrum  of  the  log  of  the  original  power  spectrum  is  called  the  power 
cepstrum.  In  the  transformation  process,  a change  is  made  from  the 
time  domain  (sec)  to  a frequency  domain  (Hz)  and  then  to  another 
temporal  domain  called  the  quefrency  domain  (sec)  which  could  be 
thought  of  as  ripples/Hz. 

If  we  consider  a<<l  then  we  can  simplify  such  that 

log  I Y(f ) | 2 - log  (1+a2)  | X(f ) 1 2-+  E lrllk+1(  2a  cos2TrfT)k 

k=l  k i+a^ 

becomes 


log|Y(f)|2  = log]x(f)  | 2+2acos2irfT 
where  the  second  term  in  this  expression  is  represented  by  the 
first  term  in  the  convergent  series.  The  power  spectrum  of  this 
expression  then  yields  the  power  cepstrum 


ypc(t)  = xpc(t)+a6(t-x) 

such  that  the  power  cepstrum  of  the  combined  wave  series  is  the 
power  cepstrum  of  the  incident  wave  series  plus  a unit  impulse 
function  or  spike  at  the  delay  time  t.  In  the  more  general  case 
when  a is  not  «1,  additional  impulses  occur  at  multiples  of  x. 

If  the  delay  time  is  sufficiently  long,  these  spikes  can 
be  observed  separated  from  the  incident  wave  power  cepstrum.  Thus 
this  procedure  can  be  followed  for  determining  the  delay  time,  x, 
between  incident  and  reflected  wave.  The  next  step  is  to  try  and 
filter  out  the  reflected  Information  found  at  x and  multiples  of  x 
to  attempt  recovery  of  the  incident  wave. 


B.  Complex  Cepstrum  and  Wave  Recovery. 


A different  approach  to  cepstral  analysis  is  seen  with  the 
investigation  of  Oppenheim  and  Schafer  into  homomorphic  deconvolution 
and  the  complex  cepstrum  {14,  7 (chap.  10)}.  Using  the  single 


reflection  model  again: 


y(t)  = x(t)+ax(t-x) 


The  Fourier  transform  of  this  is 

Y(f)  - X(f)(l+ae_2l*ifT) 

Noting  that  Y(f)  and  X(f)  are  usually  complex  quantities, 
for  this  analysis,  the  log  is  taken  directly  of  the  Fourier 
transform  of  the  time  series  and  not  of  the  power  spectrum  as  was 


done  before.  We  then  have 


log{Y  (f  ) } = log{X(f)}+log{l+ae“27rifT} 
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Using  the  log  series  expansion  expressed  previously  the 
second  term  becomes 

log(l+ae"2irifT)  = ae~12nfT-2i  e~12  (2^fT)+aie-i3(2vfT) 

2 3 

The  complete  log  expression  then  undergoes  inverse  Fourier 
transformation  to  yield  the  complex  cepstrum  for  the  combined  wave  as 

?(t)  * ft(t)+a6(t-T)^d(t-2T)+|i<S(t-3T) 

Note  that  the  complex  cepstrum  for  the  combined  wave  time  series,  £(t), 
consists  of  the  complex  cepstrum  of  the  incident  wave  time  series, 

&(t),  plus  impulse  functions  or  spikes  at  multiples  of  t which  decay 
exponentially  with  a and  also  alternate  signs.  In  all  of  the  above 
a is  considered  less  than  1.0  which  is  the  normal  characteristic  of 
the  reflection  coefficient  in  the  ocean.  If  a is  somehow  greater 
than  1.0,  a similar  examination  reveals  that  the  complex  cepstrum 
becomes  the  inverse  Fourier  transform  of  the  complex  log  of  the 
echo,  plus  impulse  functions  or  spikes  which  decay  exponentially 
with  a and  alternate  signs  in  negative  time  only  {15}.  Kemerait  and 
Childers  {15}  also  go  through  a similar  process  for  the  double 
reflection  case  which  yields  this  complex  cepstrum 
£(t)  » i(t)+a1(S(t-T1)+a26(t-T2)-^_5(t-2T1)-a1a26{t-(T1+T2)} 

(t-2t2) 

2 

where  ctj,  a2  and  T ^ , t2  represent  the  different  reflection  co- 
efficients and  delay  times  for  the  two  different  reflections  or 
echoes  from  the  same  incident  wave.  Here  ja1e“^u:'-^Tl+a2e-^ir^^T2  [ 
must  be  less  than  1.0.  This  could  be  expanded  further  to  higher 
numbers  of  multiple  echoes  from  the  same  incident  wave.  In  the 
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ocean  this  would  be  equivalent  to  recording  reflected  wave  data  from 
say  two  or  more  neighboring  breakwaters  simultaneously  combined 
with  incident  wave  data.  The  present  effort  is  limited  to  the  single 
reflection  case  only.  The  incident  ocean  wave  can  also  be  thought  of 
as  consisting  of  a summation  of  many  different  independent  waves  each 
with  its  own  single  reflection  but  more  will  be  said  on  this  later. 

In  taking  this  complex  logarithm  of  a forward  Fourier 
transform  of  a real  time  sequence,  some  important  considerations  are 
worthy  of  mention.  First,  since  with  ocean  or  water  wave  data,  only 
real  information  is  recorded,  the  Fourier  transform  of  this  real 
information  yields  complex  numbered  information.  These  complex 
numbers  are  manipula  ed  to  determine  the  phase  and  frequency  ampli- 
tude associated  with  the  sampled  wave  series,  which  becomes  a sampled 
frequency  series  after  transformation.  This  amplitude  is  found  by 
computing  the  resultant  of  each  complex  vector  and  the  phase,  <$, 
is  found  by  taking  the  inverse  tangent  of  the  imaginary  part  divided 
by  the  real  part.  In  the  x - iy  plane  this  looks  like 


Amplitude,  |x(f) j = \J (Real)2  + (Imaginary)2  ' 
Phase,  $ = tan-'*' 


The  log  of  these  amplitudes  (remembering  that  in  discrete 
analysis  each  data  point  is  acted  upon  individually)  becomes  the 
component  in  the  complex  cepstrum  computation  discussed  above.  For 
wave  recovery,  however,  the  phase  information  cannot  be  discarded. 
Additionally,  since  the  phase  thus  computed  from  complex  number 
algorithims  is  discontinuously  varying  between  -ir  and  -fir  , it  must 
be  manipulated  prior  to  taking  the  inverse  transform  as  will  be 
discussed  later. 

The  cepstrum  is  a real  valued  quefrency  (time)  series. 

As  such  it  has  no  imaginary  part.  To  achieve  this,  the  real  part  of 
the  frequency  series  undergoing  inverse  transformation  must  be  con,- 
tinuous  and  even  and  the  imaginary  part  must  be  continuous  and  odd. 


The  amplitude  of  the  frequency  series  computed  as  shown  above  as  well 
as  the  log  amplitude  will  always  yield  a continuous  and  even  function 
when  the  original  time  series  is  real.  The  phase,  however,  must  be 
made  continuous  by  "unwrapping"  and  its  linear  trend  must  be  removed 
to  avoid  undesired  shifting.  When  this  happens  the  phase  then 
becomes  continuous  and  odd.  The  log  amplitude  becomes  the  real  part 
and  the  unwrapped  phase  with  linear  trend  removed  becomes  the 
imaginary  part  of  a frequency  series  which  undergoes  inverse  Fourier 
transformation  to  yield  the  real  valued  complex  cepstrum. 

To  explain  the  phase  unwrapping  procedure  further,  Figures 
1 & 2 show  a sample  raw  phase  and  sample  unwrapped  phase.  The 
unwrapped  phase  linear  trend  is  removed  by  adjusting  data  points  so 
that  they  fall  above  and  below  a zero-mean  line  the  proportional  dis- 
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tanca  they  occur  from  a line  drawn  through  the  first  and  last  points 
of  the  unwrapped  phase  plot.  The  heavy  dotted  line  on  Figure  2 
indicates  the  line  from  which  these  distances  from  a zero-mean  are 
computed.  More  examples  of  phase  unwrapping  and  linear  trend  removal 
will  be  presented  later. 

Skinner  {11}  has  shown  that  the  power  cepstrum  can  be 
computed  from  the  complex  cepstrum  quite  easily  as 

yPC(n)  = (?(n)+?(-n))2 

That  is,  the  power  cepstrum  is  merely  4 times  the  square  of  the  even 
part  of  the  complex  cepstrum.  He  also  defines  a phase  cepstrum  that  is 
related  to  the  phase  in  the  same  way  that  the  power  cepstrum  is  related 
to  the  log  magnitude.  The  phase  cepstrum  is  four  times  the  square  of 
the  odd  part  of  the  complex  cepstrum  and  looks  like 

Y^Cn)  = (?(n)-?(-n))2 
rnC 

The  phase  cepstrum  is  presented  only  for  curiosity  as  the  power 
cepstrum  is  still  a better  indicator  of  the  delay  time,  x.  The 
phase  cepstrum  does,  however,  give  some  indication  as  to  whether 
the  phase  unwrapping  procedure  was  successful.  That  is,  if  the 
phase  cepstrum  is  very  weak  or  non  existant,  the  phase  unwrapping 
may  not  have  been  complete.  This  could  happen  if  the  phase  differed 
by  more  than  2tt  between  any  two  sample  points. 

The  procedure  for  recovery  of  the  incident  wave  first 
involves  finding  the  delay  time,  x,  by  either  the  complex  or  power 
cepstrum.  Then  the  impulse  functions  or  spikes  which  contain  the 
reflected  wave  information  found  in  the  complex  cepstrum  at  x and 


15 


multiples  of  x are  filtered  out  and  smoothed  over  by  averaging 
neighboring  points  on  either  side  of  the  spikes.  The  resulting 
filtered  complex  cepstrum  then  contains  information  pertaining  solely 
to  the  incident  wave.  By  reversing  the  process  for  computing  the 
complex  cepstrum  the  original  incident  wave  time  series  is  obtained. 
The  algorithim  for  this  process  as  obtained  from  Skinner  {16}  is 
described  by  Figure  3. 

To  demonstrate  the  wave  recovery  process,  two  models  are 
taken  from  Skinner  {11}  and  results  plotted  using  a Calcomp  Plotter. 
The  first  model  is  of  a 256-point  time  series  consisting  of  two 
64-point  long  exponential  signals  with  the  second  one  reduced  in 
amplitude  by  one  half  and  delayed  30  points.  Using  a sampling  rate 
of  .3333  sec,  this  means  the  combined  signal  is  a decaying 
exponential  pulse  21,33  sec  long  combined  with  its  echo  reduced  in 
half  and  delayed  10  sec.  Mathematically  it  looks  like: 
y(t)  = x(t)+. 5x(t-r)  where  x(t)  = te_t;  t is  actually  the  sample 
number  n multiplied  by  the  sampling  rate  It  and  T=30At=>30x. 3333  or 
x=10  sec. 

Figure  4 shows  the  raw  spectrum  of  this  combined  signal, 
and  Figure  5 shows  the  log  amplitude.  Figures  6-8  show  the  raw 
phase,  the  unwrapped  phase,  and  the  phase  with  linear  trend  removed. 
Figure  9 shows  the  real -valued  complex  cepstrum.  Figures  10  and  11 
show  the  power  and  phase  cepstrums  as  computed  from  the  complex 
cepstrum.  Note  the  large  sharp  impulse  spikes  at  multiples  of  r(10 
sec)  in  all  three  cepstral  plots  and  the  alternating  positive/ 
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negative  impulse  peaks  in  the  complex  cepstrum.  Figure  12  shows  the 
true  incident  wave  in  the  middle  with  the  recovered  incident  wave 
above  and  the  combined  wave  below.  Recovery  with  this  example  is 
almost  perfect  with  a Mean  Squared  Error  (MSE)  of  only  4.44  x 10~^ 
between  true  incident  and  recovered  incident  waves. 

The  other  signal  observed  by  Skinner  {11}  was  a 256-point 
combined  time  series  containing  two  FM  signals  pulsed  over  64  points 
of  which  one  was  reduced  by  in  amplitude  and  delayed  by  30  points 
before  being  added  to  the  incident  wave  signal.  Mathematically 
this  combined  signal  looks  like  y(t)  =»  x(t)+.5x(t-T)  where 
x(t)  “ sin ( (.4+. lit) t).  Using  a At  of  .3333  sec,  the  incident 
FM  signal  pulse  is  21.33  sec  long,  and  r=10sec.  Figures  13  through 
21  show  the  same  sequence  as  for  the  exponential  signal  but  the 
recovered  wave  similarity  to  the  incident  wave  in  Figure  21  is  much 
more  impressive  here  as  the  combined  signal  went  through  quite  a 
bit  of  unravelling.  Note  the  small  aberrations  on  the  recovered 
signal  (top  time  series  of  Figure  21)  where  not  quite  all  of  the 
reflected  wave  component  was  taken  out  probably  due  to  a bit  of 
leakage,  etc.  associated  with  the  FFT.  Mean  Square  Error  (MSE) 
between  the  recovered  incident  and  true  incident  waves  in  this 
case  was  computed  to  be  1.91  x 10“ again  indicating  almost  perfect 
recovery. 

Both  of  these  examples  were  run  under  basically  noise 
free  conditions.  Much  of  the  work  to  date  in  cepstral  analysis 
has  been  associated  with  effects  of  noise  added  to  the  signal. 
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Since  in  the  ocean  environment  wave  measurements  can  be  made  without 
noise  other  than  the  usually  negligable  amount  produced  electronically 
during  sampling  and  recording,  these  effects  will  not  be  examined  here. 
Generally  any  noise  added  to  the  signal  makes  recovery  that  much  more 
difficult. 

C.  Anticipated  Limitations  Related  to  Application  with  Ocean  Waves 

A basic  assumption  in  our  original  model  used  to  define  ceps- 
tral  analysis  is  that  the  delay  time,  t,  and  the  reflection  coefficient, 
a,  are  constant  and  not  functions  of  frequency.  With  water  waves  this 
is  not  a correct  assumption.  Both  t and  a may  vary  with  frequency  even 
in  areas  where  the  water  depth  is  constant.  When  water  depth  varies 
and  particularly  if  the  waves  under  study  are  moving  into  shallower  wa- 
ters, the  problem  is  even  more  complex  due  to  wave  celerity  varying  with 
depth  by  different  amounts  for  different  wave  frequencies.  It  is  hoped 
that  the  amount  of  differing  t's  and  a's  with  different  frequencies  will 
be  small  enough  that  this  will  not  cause  an  insurmountable  problem.  In 
this  regard  results  should  probably  be  better  for  single  peak  and  nar- 
row  banded  spectrums  than  for  double  or  multi-peaked  spectrums  with 
broader  frequency  components  as  the  different  individual  t's  will  hope- 
fully join  together  to  form  one  main  T value.  The  varying  a could  pos- 
sibly cause  the  amplitudes  of  the  Impulse  spikes  in  the  cepstrum  at 
multiples  of  t to  change  differently  than  expected  from  theory.  Both 

Although  the  literature  shows  about  equal  preference  between  the  use 
of  spectra  or  spectrums  as  the  plural  form  of  spectrum,  the  latter  form 
is  used  throughout  this  work. 
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of  these  effects  may  require  a wider  filtering  and  smoothing  of  the  re- 
flected information  spikes  in  the  complex  cepstrum  to  completely  remove 
all  reflected  wave  information  and  allow  recovery  of  the  incident  wave. 

Another  problem  that  may  be  encountered  using  the  ocean  wave 
is  the  magnitude  of  o.  Kemerait  and  Childers  {15}  found  problems  try- 
ing to  recover  the  incident  wave  as  the  reflection  coefficient  exceeds 
0.8  for  various  computer  generated  test  signals.  For  the  multiple  re- 
flection case  they  found  that  if  the  reflected  signal  amplitudes  are 
less  than  the  basic  wave,  recovery  is  possible.  In  the  general  wave 
case,  if  the  reflected  signal  amplitudes  are  greater  than  the  incident 
wave,  but  not  equal,  the  largest  reflected  signal  can  be  recovered.  If 
two  or  more  reflected  signals  have  larger  amplitude  than  that  of  the 
incident  wave,  recovery  is  not  possible.  Since  ocean  structures  fre- 
quently may  have  reflection  coefficients  approaching  1.0,  this  real 
limitation  will  be  examined  for  at  least  the  single  reflection  case. 

Cepstral  analysis  to  date  has  been  involved  with  pulsed  and 
burst-type  signals  of  relatively  short  duration  similar  to  the  recent- 
ly discussed  exponential  and  FM  signal  examples.  A problem  may  arise 
from  using  combined  incident/reflected  wave  data  collected  over  the 
total  sample  length.  Since  the  time  series  is  in  fact  digitally  sam- 
pled it  is  felt  this  should  not  be  an  insurmountable  problem.  It  may 
mean  zeroes  will  have  to  be  added  to  the  end  of  the  sampled  data 
which  increases  the  frequency  sampling  rate  to  get  results  similar  to 
that  of  the  pulsed  signal.  More  will  be  said  about  this  in  Chapter  V. 
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CHAPTER  IV 

DESCRIPTION  OF  TEST  MODELS 

A.  Summation  of  Finite  Number  of  Sinusoids  with  Random  Phase. 

The  first  test  model  for  ocean  simulated  waves  maintains 
both  a,  the  reflection  coefficient,  and  x,  the  delay  time,  as 
constants  non-variant  with  respect  to  frequency.  The  incident  wave 
model  used  in  this  test  is  that  derived  by  Borgman  {17}  which  simulates 
an  ocean  waveform  based  on  an  assumed  spectral  curve.  This  model  can 
be  described  as  the  summation  of  a finite  number  of  sinusoid  components 
having  random  phase  shift  angles.  Simplified  mathematically  it  looks 
like: 

x(t)  =»  I — c r2*  cos  (2irf  At+A  ) 

m-1  V M 31  m 

where  x(t)  represents  the  simulated  ocean  wave  surface  profile.  As  M, 
which  represents  the  number  of  sinusoid  components,  gets  larger,  the 
more  realistic  the  simulated  wave  becomes.  The  tradeoff  for  using 
larger  and  larger  values  of  M is  increased  computer  cost  for  "building" 
the  simulated  wave.  For  the  testing  in  this  work,  an  M value  of  40  is 
used  for  a 256-point  time  series.  The  amplitude  of  each  component  is 
found  by  dividing  the  total  assumed  energy  represented  as  the  var- 
iance, a2,  of  the  assumed  spectrum  (where  a2  is  the  area  under  that 
spectral  curve)  by  the  number  of  sinusoids  being  summed  and  taking  the 
square  root  of  that  result  (remembering  that  the  square  root  of 
energy  is  directly  proportional  to  wave  height).  Therefore,  the  total 
energy  in  the  simulated  wave  will  equal  the  energy  represented  by  the 
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initially  assumed  spectral  curve.  When  the  assumed  spectral  curve  is 
divided  into  equal  area  increments,  fm  values  are  chosen  as  the  mid- 
frequency of  each  incremental  frequency  band.  At  is  the  time  increment 
between  samples,  and  (J>m  is  a uniformly  distributed  random  variable 
between  zero  and  2ir  which  represents  the  random  phase  shift  or 
starting  point  of  each  sinusoid  component  being  summed. 


For  example,  with  an  M value  of  5 the  assumed  spectrum  would 
be  divided  as  shown  below  to  obtain  the  equal  energy  value  for  comput- 
ing amplitudes  and  for  determining  f values  or  the  mid-point  fre- 
quencies . 


(not  to 
scale) 


frequency 


The  assumed  spectrum  used  in  creating  this  wave  form  is 
that  proposed  by  Pierson  and  Moskovitz  {18}  and  modified  for  limited 
fetches  by  Mitsuyasu  {19},  Barnett  {20}  and  Sylvester  and 


12  3 4 5 


21 


is  the  wind  velocity,  f the  frequency,  and  g the  gravitational 
constant.  However,  as  derived  by  Christensen  and  described  by 
Morden  {2},  for  limited  fetches  a and  3 become 


a 


5(2ir)Vf£  - 3.49(2^r)4S(fp)fp5 
g2  g2 


3 


5(2TrUf£)4 
^ g 


where  f is  the  frequency  of  the  spectral  peak  and  S(fp)  is  the 
amplitude  of  the  spectral  peak.  The  variance  a2  is  then  found  from 

a2  = 2/°°S(f  )df  = a^, 

° 4Bg2 

A model  wave  then  can  be  created  based  on  the  variables  U, 
fp,  S(fp),  and  a2.  Assuming  U and  fp  and  either  a2  or  S(fp)  will 
completely  define  the  P-M  spectrum. 

Therefore,  the  simulated  incident  ocean  wave  is  computed 
by  first  assuming  values  of  S(fp)  and  f to  compute  the  P-M  spectral 
model  and  a value  for  a2.  This  spectral  curve  is  divided  into  equal 
areas  and  mid-point  frequencies  for  each  equal  area  increment  are 
computed.  Using  a suitable  random  number  generator  to  compute  random 
phase  angles,  each  sinusoid  is  computed  and  summed  to  the  previously 
computed  sinusoid (s)  as  described  previously.  For  simplicity,  this 
simulated  incident  ocean  wave  will  be  called  the  P-M  generated  wave. 

Taking  the  P-M  generated  wave  and  scaling  it  by  a 
reflection  coefficient  a produces  the  simulated  reflected  wave. 
Delaying  this  reflected  wave  by  x and  adding  the  two  together  yields 
the  initial  P-M  generated  combined  wave  test  model.  Note  x is 
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constant  for  the  whole  wave  train  in  this  initial  model.  In  Chapter  V, 
further  description  of  model  manipulation  is  described  related  to  test- 
ing procedure. 

3.  Varying  Delay  Time  with  Respect  to  Frequency 


In  water  waves,  the  delay  time  x varies  with  respect  to  the 

wave  celerity  which  varies  as  a function  of  wave  length  which  in  turn 

varies  with  frequency,  so  it  follows  that  x does  in  fact  vary  with 

frequency.  For  linear  deep  water  wave  theory  this  variation  can  be 

xf 


expressed  in  fps  units  by  x 


2.56 


where  f is  the  frequency  of  the 


wave  component  and  x is  the  distance  from  the  wave  measurement  device 


to  the  reflecting  surface.  This  relationship  is  derived  from 

distance  traveled 

celerity  = ; " = 

time  of  travel 

where  T is  the  wave  period  and  g is  the  gravitational  constant. 


— £•  T for  deep-water  wave  theory 
2tt 


2v  i 

Therefore  celerity  * — = 5.12  — since  distance  x is  traveled  twice 

T f 


over  the  delay  time,  x,  so  that  x 


2xf 


xf 


5.12  2.56 

Using  the  P-M  generated  wave  model  previously  described,  the 
incident  wave  is  generated  in  the  same  manner  as  before.  However,  the 
reflected  wave  components  are  delayed  by  their  own  x values  for  their 
particular  frequency  before  being  summed  together,  scaled  by  a and 
combined  with  the  incident  wave.  This  is  then  a model  where  each  sep- 
arate reflected  frequency  component  has  its  own  delay  time  x contri- 
buting to  the  overall  resultant  x of  the  reflected  wave  time  series. 
Results  from  using  this  model  are  discussed  in  the  next  chapter. 
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C.  Varying  Reflection  Coefficient  with  Respect  to  Frequency. 

To  carry  the  P-M  generated  wave  model  one  step  closer  to 
the  real  world  reflected  water  wave,  an  attempt  is  made  at  varying 
the  reflection  coefficient  as  a function  of  frequency.  This  first- 
try  variation  is  performed  by  merely  making  a a linear  function  vary- 

. 2 

ing  from  .4  to  .6,  .5  to  .7,  or  .7  to  -.9  incremented  by  for  each 
sinusoidal  component  of  the  reflected  wave  where  M is  the  total  number 
of  components.  The  results  from  using  this  test  model  are  also  dis- 
cussed in  Chapter  V. 

In  reality,  the  reflection  coefficient,  ct,  probably  varies 
as  some  non-linear  function  of  frequency  such  that  a would  be  higher 
for  higher  frequency  waves  than  for  lower  frequency  waves  depending  on 
the  structure.  That  is,  the  lower  frequency  wave  components  are  more 
likely  to  be  absorbed  or  transmitted  by  most  structures  than  are  the 
higher  frequency  wave  components  making  up  the  ocean  wave.  However, 
the  simplified  model  described  above  and  in  Chapter  V is  sufficient 
for  an  initial  test  of  effects  from  a varying  as  a function  of 
frequency. 


CHAPTER  V 


APPLICATION  OF  CEPSTRAL  ANALYSIS  TO  TEST  MODELS 

As  mentioned  in  Chapter  III,  cepstral  analysis  techniques 
are  usually  applied  to  "pulse"  type  signals  where  no  signal  informa- 
tion is  contained  between  the  incident  wave  and  the  reflected  wave. 
Figure  22  shows  a sample  of  incident  and  reflected  pulse  signals  which 
would  be  added  together  to  form  such  a combined  incident/reflected  sig- 
nal. Note  the  lack  of  information  over  time,  x,  at  the  end  of  the 
incident  wave  and  over  time,  t,  at  the  beginning  of  the  reflected  wave. 
Both  incident  and  reflected  waves  are  of  the  same  length.  The  expon- 
ential and  FM  signals  discussed  in  Chapter  III  are  examples  of  pulse 
waves . 

When  the  incident  wave  extends  the  full  length  of  the  time 
series,  the  incident  and  reflected  waves  appear  as  shown  in  Figure  23. 
This  combined  wave  is  called  a "pulse-shifted"  wave  as  no  infor- 
mation exists  over  time,  t,  in  the  beginning  of  the  reflected  wave. 

Note  the  reflected  wave  is  now  not  as  long  as  the  incident  wave.  This 
creates  a situation  known  as  "echo  truncation"  where  the  incident  wave 
does  not  terminate  x seconds  before  the  end  of  the  reflected  wave  as 
in  the  pulse  signal  described  above. 

For  a more  realistic  look  at  ocean  wave  data  sampling,  a 
"forward-shift"  combined  wave  is  created  from  incident  and  reflected 
waves  following  the  pattern  shown  in  Figure  24.  The  name  "forward- 
shift"  comes  from  the  fact  that  a new  incident  wave  is  made  from  the 
previously  created  incident  wave  by  extending  the  wave  train  and 


figuratively  shifting  it  forward.  The  important  item  here  is  that 
signal  information  exists  from  both  the  incident  and  reflected  waves 
throughout  the  sampling  period. 

To  test  the  first  model  described  in  Chapter  IV  where  neither 
a nor  x are  allowed  to  vary  with  frequency  the  combined  wave  is  formed 
first  using  the  "pulse-shift"  combination  and  then  using  the  "forward- 
shift"  combination.  In  both  cases,  the  incident  and  reflected  waves 
are  created  and  stored  and  each  data  point  in  the  reflected  wave  is 
multiplied  by  a reflection  coefficient  before  the  two  waves  are  added 
together  to  form  the  combined  wave. 

The  P-M  Spectral  model  used  to  create  the  P-M  Generated  Wave 
is  shown  as  Figure  25.  This  spectrum  is  selected  as  most  of  the  spec- 
tral energy  is  concentrated  around  a narrow  frequency  band.  With  later 
trials  this  should  help  produce  one  x value  for  the  total  combined  wave 
when  x is  allowed  to  vary  with  frequency. 

Using  a rectangular  window  (no  weighting  of  input  data)  as 
recommended  by  Skinner  {11}  for  the  pulse-shifted  wave,  results  are 
somewhat  discouraging.  For  a 256  data  point  time  series.  At". 3333  sec, 
x=10  sec  (30  data  points)  and  a=.8,  the  Mean  Squared  Error  (MSE)  between 
the  recovered  incident  wave  and  the  true  incident  wave  is  .122  with  the 
MSE  between  the  true  incident  and  combined  wave  equal  to  .285.  The  MSE 
(equal  to  the  mean  of  the  difference  squared  comparing  parallel  data 
points  of  two  "side-by-side"  time  series)  between  recovered  and  true 
incident  waves  is  an  indicator  of  recovery  results  and  the  MSE  between 
true  incident  and  combined  waves  is  an  indicator  of  what  the  recovery 
process  starts  with.  Comparing  these  two  indicators,  the  MSE  relation- 
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ships  show  that  some  recovery  did  occur  but  that  the  recovered  incident 
wave  and  the  true  incident  wave  are  still  not  very  much  alike. 

Oppenheim  and  Shafer  {7}  state  that  using  an  exponential  win- 
dow (exponential  weighting  of  input  data)  enhances  the  identification 
of  delay  time,  t,  and  improves  the  overall  cepstral  analysis  recovery 
process.  Skinner  {11}  found  this  not  to  be  the  case  for  the  exponen- 
tial and  FM  signals  discussed  in  Chapter  III.  Using  an  exponential 
window  on  the  pulse-shifted  combined  wave  where  each  data  point  (n=l  to 
n=N)  is  "weighted"  or  multiplied  by  (.99)n-^  before  the  cepstral  decom- 
position process  and  dividing  the  resulting  recovered  wave  by  (.99)°^ 
to  "unweight"  the  output  lowers  the  MSE  between  original  and  recovered 
incident  waves  to  .044.  This  indicates  quite  satisfactory  recovery  re- 
sults as  can  be  observed  from  Figure  34. 

It  should  be  pointed  out  here  that  although  "one-sided"  plots 
(representing  only  the  positive  half  of  the  spectrum)  of  the  raw  spec- 
trum and  log  amplitude  are  shown  for  the  exponential  and  FM  signals, 
all  subsequent  plotting  of  raw  spectrum  and  log  amplitude  before  the 
cepstrum  is  computed  are  shown  as  "two-sided"  plots  (noting  that  the 
negative  half  of  the  spectrum  is  reproduced  to  the  right  of  the  posi- 
tive half  since  transformation  was  performed  via  the  DFT) . This  is 
done  to  show  that  in  fact  the  log  amplitude  is  indeed  continuous  and 
even  and  the  unwrapped  phase  with  linear  trend  removed  is  continuous 
and  odd  which  together  produce  a real  valued  cepstrum.  Smoothed  spec- 
tral curves  to  be  seen  later  will  still  be  shown  as  one-sided  plots. 

To  discuss  the  positive  results  from  exponential  weighting 
further,  it  is  seen  that  it  helps  reduce  echo  truncation  error  by  plac- 
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ing  greater  emphasis  on  the  beginning  of  the  input  combined  wave  time 
series  than  on  the  end.  Different  exponential  "base"  values  from  .900 
to  .999  were  tried  for  different  input  conditions  and  a base  value  of 
.99  always  seemed  to  produce  the  best  results. 

Figures  26  through  34  show  the  plots  of  the  time  series 
manipulation  through  the  spectral  and  cepstral  processes  at  the  same 
stages  as  listed  for  the  exponential  and  FM  signals.  An  important  ob- 
servation when  seeing  these  plots  is  that  this  more  irregular  wave  form 
produces  a more  irregular  cepstrum  but  the  reflected  information  im- 
pulse spikes  are  stil,i  very  obvious  in  all  three  cepstrums.  When  the 
spikes  located  at  multiples  of  t are  filtered  out  of  the  complex  ceps- 
trum, wave  recovery  is  quite  good  as  is  seen  in  Figure  34  when  compar- 
ing the  original  incident  time  series  located  in  the  middle  of  this 
figure  to  the  recovered  incident  time  series  located  at  the  top. 

A better  indicator  of  how  good  the  results  are  using  this 
model  is  seen  in  Figures  35  through  38.  These  figures  represent  the 
smoothed  spectrums  (every  8 points  averaged  together  for  these  parti- 
cular plots)  of  the  combined  wave,  the  original  incident  wave,  the  re- 
covered incident  wave,  and  the  reflected  wave,  respectively.  The  spec- 
trums of  the  recovered  and  original  incident  waves  are  almost  identical 
indicating  good  wave  recovery.  The  reflected  wave  spectrum  is  obtained 
from  the  recovered  reflected  wave  found  by  subtracting  the  recovered 
incident  wave  data  points  from  the  combined  wave  data  points.  Using 

the  formula  a?  = a 2a^  which  relates  the  variance,  or  area  under  the 
K I 


spectrum  of  the  reflected  wave  to  the  variance  or  area  under  the  spec- 
trum of  the  incident  wave  times  the  reflection  coefficient  squared,  an 
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attempt  is  made  at  determining  the  reflection  coefficient  used  to  cre- 
ate the  combined  wave.  The  resulting  computation  produces  a reflection 
coefficient  equal  to  .785,  which  is  very  close  to  the  actual  .8  value 
used. 


The  results  of  this  first  test  of  a P-M  generated  wave  shape 
therefore  are  very  encouraging.  The  next  tests  using  the  "forward- 
shift"  model,  however,  are  not  so  gratifying.  Input  parameters  are  the 
same  as  with  the  pulse-shift  model.  Figures  39  through  47  show  the 
same  plot  sequence  as  before.  The  impulse  spikes  are  not  as  readily 
noticable  in  the  three  cepstrums  now  and  recovery  attempts  as  seen  from 
the  time  series  comparisons  in  Figure  47  between  the  original  incident 
wave  (middle)  and  the  recovered  incident  wave  (top)  show  results  poorer 
than  with  the  pulse-shift  test.  The  smoothed  (averaged  spectral  value 
of  every  10  data  points)  spectral  plots  of  results  shown  in  Figures  48 
through  51  also  show  recovery  as  less  than  satisfactory.  The  recovered 
wave  spectrum  shows  its  energy  level  to  be  somewhere  between  that  of 
the  true  incident  and  combined  waves  instead  of  equal  to  that  of  the 
true  incident  wave.  The  reflected  wave  spectrum  is  obviously  too  small 
such  that  the  computed  reflection  coefficient  is  much  lower  than  the 
value  of  .8  used  in  creating  the  combined  wave.  An  encouraging  obser- 
vation in  this,  however,  is  that  some  energy  is  in  fact  removed  from 
the  combined  wave  so  that  the  process  at  least  attempted  to  perform 
some  decomposition. 

Following  Skinner's  {11}  recommendation,  a first  step  or  try 
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at  improving  results  is  to  add  N zeroes  to  the  end  of  the  input  data 
(adding  256  zeroes  for  256  data  points)  before  performing  the  spectral 


and  cepstral  computations.  This  effectively  doubles  the  sampling  rate 
in  the  frequency  domain  and  hence  the  quefrency  domain  without  changing 
the  time  series  information  whatsoever.  When  zeroes  are  added  to  the 
end  of  the  pulse-shifted  combined  wave  from  before,  no  improvement  is 
noted,  probably  due  to  the  fact  that  results  are  already  very  good. 
Adding  N zeroes  to  the  end  of  the  input  data  of  the  forward-shifted 
combined  wave  before  performing  decomposition  yields  the  smoothed  re- 
covered incident  and  reflected  wave  spectral  plots  shown  in  Figures 
52  and  53.  It  is  seen  that  this  technique  helps  produce  a recovered 
incident  wave  spectrum  that  is  slightly  closer  to  the  true  incident 
wave  spectrum,  but  results  are  still  unsatisfactory  as  evidenced  by  the 
again  small  reflected  wave  spectrum  and  the  obvious  resulting  incorrect 
calculation  of  the  reflection  coefficient . 

Manipulation  of  data  sampling  methods  is  also  employed  in 
searching  for  better  results.  These  methods  include  trying  different 
time  sampling  rates,  trying  different  r values  and  increasing  the  total 
sample  length  to  1024  data  points  from  the  previous  length  of  256  data 
points.  All  attempts  at  data  sampling  manipulation  seem  to  be  ineffec- 
tual for  improving  results.  When  the  reflection  coefficient  is  reduced 
to  .5,  results  improve  only  very  slightly.  Attempts  at  improving  re- 
covery using  a rectangular  window  and  an  exponential  window  with  dif- 
fering base  values  above  and  below  .99  show  best  results  are  still 
achieved  with  exponential  weighting  and  un-weighting  using  the  (.99)n-^- 
factor  described  previously. 

Therefore,  a problem  exists  using  the  P-M  Generated  wave 
forward-shift  model  even  when  t and  a are  constant  and  independent  of 
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frequency.  This  problem  appears  to  be  centered  around  the  fact  that 
information  is  now  contained  in  the  beginning  of  the  reflected  wave. 
This  information  looks  more  like  noise  when  comparing  the  incident  and 
reflected  waves’  time  series  because  the  actual  image  reflection  from 
the  incident  wave  doesn't  start  until  t seconds  later. 

When  the  second  model  described  in  Chapter  IV  is  employed, 
where  x is  allowed  to  vary  with  frequency  for  each  sinusoid  making  up 
the  P-M  generated  combined  wave,  even  less  incident  wave  recovery  is 
possible.  The  different  attempts  at  improving  results  follow  the 
same  sequence  as  those  tried  with  the  forward  shift  wave.  Figure  54 
shows  a typical  complex  cepstrum  for  this  model  where  At=.3333,  x=10 
sec  and  ct=.5.  It  is  readily  obvious  that  identifying  a x,  even  when 
the  input  x is  known,  is  very  difficult  in  this  situation.  Filtering 
the  impulse  spikes  which  should  be  present  at  multiples  of  x yields 
the  filtered  complex  cepstrum  shown  in  Figure  55.  Figures  56  and  57 
demonstrate  the  same  technique  showing  a complex  cepstrum  and  a 
filtered  complex  cepstrum  for  At=.7333.  Again  it  is  quite  difficult 
to  identify  ax. 

Since  it  appears  initially  that  the  x peak  in  Figure  54  may 
be  wider  than  just  one  point,  five  different  filter  variations  are  at- 
tempted filtering  out  more  than  one  point.  Some  examples  are  shown  in 
Figures  58  and  59  of  the  resulting  filtered  complex  cepstrums.  None 
of  these  different  filter  methods  help  and  in  fact  all  of  them  make 
the  MSE  even  greater  between  the  original  and  recovered  incident  waves. 
Some  attempts  even  made  the  MSE  between  actual  incident  and  recovered 
incident  waves  higher  than  between  the  incident  and  combined  waves. 
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Therefore,  it  appears  the  impulse  spike  definition  for  location  of  t is 
more  accurate  than  a wider  peak  definition  of  t even  in  the  model  where 
r is  allowed  to  vary  with  frequency. 

When  employing  the  third  model  described  in  Chapter  IV,  where 
reflection  coefficient  is  also  allowed  to  vary  with  frequency,  very 
little  difference  is  observed  in  the  appearance  of  the  complex  cepstrum 
using  the  same  input  parameters  as  for  the  constant  reflection  coef- 
ficient model  above.  Wave  recovery  attempts  for  this  model  are  not  sa- 
tisfactory. Again,  no  improvement  is  obtained  by  using  different  At's, 
different  t's,  or  different  reflection  coefficient  ranges,  including  a 
varying  from  .7  to  .9,  .5  to  .7,  and  .4  to  .6.  Example  complex  ceps- 
trums  for  this  model  are  shown  in  Figures  60  and  61.  The  difficulty  in 
identifying  a value  for  t is  also  readily  apparent  for  this  model. 

The  problems  with  the  latter  two  models  cannot  be  readily 
identified  until  satisfactory  results  can  first  be  obtained  from  the 
forward  shift  model  where  a and  t are  non-variant  with  frequency. 

This  requires  some  way  of  satisfactorily  dealing  with  the  "noise"  or 
continuous  wave  section  existing  in  front  of  the  reflected  wave  con- 
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CHAPTER  VI 

APPLICATION  OF  CEPSTRAL  ANALYSIS  TO  ACTUAL  OCEAN  DATA 

Some  time  series  data  records  recorded  near  a breakwater  in 
Sitka,  Alaska,  in  1973  are  run  through  this  algorithm  to  determine  if  a 
value  for  t can  be  identified.  Figure  62  shows  the  breakwater  plan  and 
the  location  of  instruments  on  this  breakwater.  The  data  points  util- 
ized are  taken  from  the  combined  wave  recorded  at  the  Northernmost 
gage  when  wave  action  was  predominately  perpendicular  to  the  short  leg 
of  the  breakwater.  This  site  has  the  advantage  that  any  recovered  in- 
cident wave  obtained  through  the  decomposition  process  can  be  compared 
for  accuracy  with  the  true  incident  wave  recorded  from  the  wave  staff 
off  the  corner  of  the  breakwater. 

These  data  records  had  a major  disadvantage,  however,  as 
the  breakwater  is  only  about  two  feet  deep  allowing  some  of  the  low 
frequency  waves  to  transmit  through  the  breakwater  causing  the  reflec- 
ted wave  to  be  different  in  shape  and  frequency  composition  from  the 
incident  wave.  Also,  the  frequency  spectrum  for  this  location  is 
"double-peaked,"  that  is,  a low  frequency  swell  condition  exists  as 
well  as  the  waveform  produced  by  local  winds  over  a local  fetch. 

Looking  at  cepstrums  generated  from  this  data  is  informative, 
however,  in  that  it  reveals  some  of  the  other  problems  associated  with 
using  this  process  on  real  world  data.  Three  records  are  selected  and 
tested  which  are  2048  data  points  long  recorded  at  a At  of  .44  sec. 
producing  15  minute  long  data  records.  An  example  of  the  complex 
cepstrum  thus  obtained  using  no  weighting  of  the  input  data 
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(rectangular  window)  is  shown  in  Figure  63..  It  is  interesting 
to  observe  the  unwrapped  phase  and  unwrapped  phase  with  linear  trend 
removed  for  this  data  which  are  shown  in  Figures  64  and  65.  The 
simularity  between  the  two  raises  suspicion  that  phase  unwrapping 
may  not  have  been  complete.  A remedy  for  this  problem  is  to  add 
zeroes  to  the  end  of  the  original  input  data.  This  is  tried  using 
1024  data  points  from  the  same  record  and  adding  1024  zeroes  before 
computing  the  spectrum  and  complex  cepstrum. 

Figures  66,  67,  68,  and  69  show  the  raw  spectrum,  log 
amplitude,  unwrapped  phase,  and  unwrapped  phase  with  linear  trend 
removed  for  this  new  input.  The  linearity  seen  in  the  unwrapped  phase 
plot  shows  phase  unwrapping  is  indeed  now  more  complete. 

Figures  70  and  71  show  the  complex  cepstrum  and  power 
cepstrums,  respectively.  Difficulty  in  obtaining  a value  for  t from 
these  plots  is  readily  apparent.  Perhaps  some  help  could  ie  obtained 
by  plotting  the  cepstrum  over  a longer  quefrency  axis  for  better  plot 
visualization.  Using  deep  water  wave  theory,  the  physical  location 
of  the  wave  staff  and  predominant  frequency  from  the  wave  spectrum, 
x is  estimated  to  be  about  18-20  seconds.  This  would  probably  be 
very  hard  to  identify  in  a 15  minute  record  even  if  an  impulse  spike 
were  located  at  this  point. 

The  most  important  reason  for  failure  to  identify  a i for 
this  situation  is  more  probably  due  to  the  fact  that  the  reflected 
wave  is  altered  in  shape  and  frequency  and  is  not  in  fact  a true 
replica  of  the  original  incident  wave.  It  seems  Cepstral  Analysis 
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decomposition  should  work  even  if  the  reflected  wave  shape  is 
different  from  the  incident  wave  shape  as  long  as  the  frequency 
information  is  not  altered  in  the  reflection  process.  Since  the 
floating  breakwater  used  here  does  in  fact  filter  certain  low 
frequency  components  out  of  the  reflected  wave  by  allowing  them  to 
pass  through  the  breakwater,  it  is  not  surprising  that  an  impulse 
spike  is  not  observable  in  either  the  complex  or  power  cepstrums  at 
an  identifiable  location  for  x. 

The  predominant  spike  at  the  beginning  of  these  cepstrums 
is  left  over  from  the  incident  wave  complex  cepstrum  which  always 
produces  an  extremely  large  spike  at  zero  quefrency  in  the  cepstral 
plots.  In  fact,  all  of  the  cepstral  plots  seen  in  this  thesis  have 
their  first  5 data  points  set  to  zero  amplitude  on  either  side  of 
zero  quefrency  (more  properly  in  the  digital  sense,  points 
n - 1,  2,  3,  4,  5,  N,  N-l,  N-2,  N-3,  and  N-4  were  set  equal  to  zero) 
to  allow  easier  identification  of  the  smaller  amplitude  spikes 
containing  the  reflected  wave  information.  Otherwise,  plotting 
would  be  based  on  the  large  spike  at  zero  quefrency  and  the  smaller 
reflected  information  impulse  spikes  would  be  even  more  difficult 
to  identify. 

Even  taking  this  into  consideration,  an  actual  x value 
is  practically  impossible  to  identify  from  these  cepstrums.  The 
power  cepstrum  shown  in  Figure  71  offers  many  possible  impulse  spikes 
but  at  times  which  are  totally  unrealistic  with  the  ocean  wave  data 
interacting  with  this  breakwater.  Obviously,  better  data  and  more 
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CHAPTER  VII 

CONCLUSIONS  AND  RECOMMENDATIONS 

Cepstral  analysis  decomposition  techniques  work  extremely 
well  with  the  pulse-snift  P-M  generated  combined  wave  model  and  not 
so  well  with  the  forward-shift  P-M  generated  combined  wave  model  as 
described  in  Chapters  TV  and  V.  Problems  exist  with  identification 
of  delay  time  as  well  as  with  decomposition  of  the  combined  wave  when 
cepstral  analysis  techniques  are  applied  to  P-M  generated  wave 
models  which  vary  t and/or  the  reflection  coefficient,  a,  as  a 
function  of  frequency.  More  research  is  required  before  this 
decomposition  technique  can  be  successfully  applied  to  real-world 
ocean  waves.  It  may  never  work  for  decomposition  of  combined 
waves  for  which  the  reflected  wave  frequency  information  is 
altered  by  the  structure  and  no  longer  is  the  same  as  the  incident 
wave,  such  as  occurs  with  shorelines  and  many  breakwaters. 

Exponential  windowing  or  weighting  does  help  when  using 
these  decomposition  techniques  with  P-M  generated  wave  models. 

Adding  zeroes  to  the  end  of  the  input  data  prior  to  application  of 
cepstral  analysis  decomposition  does  improve  phase  unwrapping  but 
helps  only  slightly  to  none  at  all  with  incident  wave  recovery. 

Of  the  three  cepstrum  types  examined,  the  power  cepstrum 
is  the  best  indicator  of  delay  time  between  incident  and  reflected 
waves.  The  phase  cepstrum  does  not  seem  to  offer  much  information 
other  than  reinforcement  or  possible  validation  of  the  power 
cepstrum  information. 
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It  is  recommended  that  this  process  be  tested  using  real- 
world  ocean  data  for  which  the  reflected  wave  frequency  information 
is  identical  to  that  of  the  incident  wave  such  that  all  frequency 
components  are  reflected  by  the  structure.  Further  experimentation 
should  begin  using  real-world  data  with  a "single-peaked"  smoothed 
spectrum  before  going  to  data  with  a "two-peaked"  or  "multi-peaked" 
spectrum.  When  using  data  with  more  than  one  spectral  peak, 
experimentation  could  be  made  into  the  effect  of  filtering  out  one 
or  more  peak(s)  before  going  through  the  decomposition  process.  Of 
particular  interest  would  be  results  obtained  after  filtering  out 
the  lower  frequencies  from  the  combined  wave  that  would  be 
transmitted  through  a structure  like  the  floating  breakwater 
investigated  in  Chapter  VI. 

Since  the  complex  cepstrum  consists  of  impulse  spikes 
at  multiples  of  t which  decay  exponentially  as  a function  of  the 
reflection  coefficient,  investigation  could  be  made  into  possible 
determination  of  reflection  coefficients  directly  from  the 
amplitudes  of  these  peaks.  This  would  save  the  inverse  process 
presently  necessary  by  not  requiring  actual  determination  of  the 
incident  wave  to  obtain  the  reflected  wave  before  computing  a 
reflection  coefficient. 

Research  should  be  conducted  towards  improving  results 
with  the  forward-shift.  One  possible  alternative  is  to  try 
reiterative  techniques.  Once  a value  for  delay  time  is  determined, 
the  recovered  wave  could  possibly  be  improved  by  putting  it  through 
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the  decomposition  process  again.  Another  possibility  to  improve 
results  is  to  somehow  weight  or  alter  the  data  points  at  the 
beginning  of  the  combined  wave  time  series  so  that  the  leading 
portion  of  the  reflected  wave  will  not  have  such  a degrading  effect 
on  results. 

It  appears  from  attempts  at  adjusting  the  filtering  of 
the  impulse  spikes  from  the  complex  cepstrum  that  the  spikes  do 
not  join  together  to  form  a wider  peak  when  r varies  with  frequency. 
Perhaps  what  actually  happens  is  impulse  spikes  occur  separately 
at  the  different  t values  produced  by  different  frequency 
components  making  up  the  wave.  This  possibility  could  be 
investigated  further. 

Other  windowing  or  weighting  of  input  data  prior  to  the 
decomposition  process  has  generally  not  met  with  success  when 
using  cepstral  analysis  techniques.  However,  this  could  be 
investigated  further  for  the  ocean  wave  by  experimenting  with  other 
data  windows. 

Generally,  cepstral  analysis  decomposition  methods  do  not 
appear  at  this  time  to  be  the  final  answer  in  recovery  of  the 
incident  wave  from  combined  incident/reflected  ocean  wave  data. 

First,  the  r value  itself  must  be  made  reliably  discernable  in  the 
cepstrum  and  then  actual  recovery  must  be  made  more  efficient. 

However,  with  further  research,  it  may  yet  prove  to  be  a valuable  tool 
in  helping  to  determine  delay  times,  reflection  coefficients,  and 
other  parameters  associated  with  wave  decomposition  when  determining 
true  wave  forces  acting  on  and/or  from  a structure. 
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Figure  1.  Sample  Raw  Phase  Plot 


Sample  Unwrapped  Phase  Plot 


Figure  3.  Wave  Recovery  Algorithm 
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Recovered  Incident,  True  Incident  and  Combined  Wave  Time  Series  for  FM  Signal 
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Figure  25.  Original  P-M  Spectral  Model 
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Figure  27.  Log  Amplitude  of  Pulse-Shifted  Combined  P-H  Generated  Wave 
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Recovered  Incident,  True  Incident  and  Combined  Wave  Time  Series  for  Pulse-Shifted 
P-M  Generated  Wave 
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Figure  38. 
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Smoothed  Spectrum  of  Reflected  P-M  Generated  Wave  from 
Pulse-Shifted  Combined  P-M  Generated  Wave 
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Figure  44.  Complex  Cepstrum  for  Forward-Shifted  Combined  P-M  Generated 
Wave 
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Figure  47.  Recovered  Incident,  True  Incident  and  Combined  Wave  Time  Series  for  Forward 
Shifted  P-M  Generated  Wave 
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49.  Smoothed  Spectrum  of  Incident  P-M  Generated  Wave 
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>1.  Smoothed  Spectrum  of  Reflected  P-M  Generated  Wave  from 
Forward-Shifted  Combined  P-M  Generated  Wave 
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Figure  55.  Filtered  Complex  Cepsti 
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Figure  62.  Sitka,  Alaska,  Breakwater  Plan  and  Instrumentation  Layout 
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Data  Record  of  Sitka  Breakwater  Combined  Wave  with  1024 
Zeroes  Added  to  End 
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Figure  70. 
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Complex  Cepstrum  of  1024-Point  Data  Record  of  Sitka 
Breakwater  Combined  Wave  with  1024  Zeroes  Added  to  End 


